function [ ] = plotFitScatterEmissionsAll( plant, ctrM, scenarios_toplot, ...
    filestub, titlestub, inlogs, optfig )
% Plot total costs under the emissions market and under the command and
%   control scenario
%
%   INPUTS:
%     model - Model object with simulated data
%     ctrM  - Counterfactual object with data simulated from estimated
%             parameters
%     scenarios_toplot - Table of scenarios to plot
%     filestub - Root of filename in which to print figure
%     titlestub - Root of title to give the figure
%     inlogs - Whether to plot emissions histogram in logs
%     optfig - Figure options

% Pick out only the scenarios passed as arguments to plot
scenarios_toplot = renamevars(scenarios_toplot,'cap','id_level');
ctrTable = ctrM.output_itk;
ctrTable = innerjoin(ctrTable,scenarios_toplot,...
      'Keys',{'id_period','id_level'});
tabulate(ctrTable.id_period);

% Get axis limits across all periods
Eitall = [ plant.Eit; ctrTable.Eit ];
Elim = prctile(log(Eitall),[0.1 99.9]);
Elim(1) = floor(Elim(1));
Elim(2) = ceil(Elim(2));

% Plot scatter of emissions for each period
for i = 1:size(scenarios_toplot,1)
    
    % Get parameters describing the scenario: period, cap
    t   = scenarios_toplot.id_period(i);
    cap = scenarios_toplot.id_level(i);
    
    % In data (model)
    keepd = (plant.id_period==t & plant.treatment==1);
    plant_it = plant( keepd, {'id_plant','Eit'} );
    plant_it = renamevars(plant_it,'Eit','Eit_data');
    
    % In simulation (counterfactual)
    keeps = (ctrTable.id_period == t & ctrTable.id_level == cap);
    model_it = ctrTable( keeps, {'id_plant','Eit'} );
    model_it = renamevars(model_it,'Eit','Eit_model');
    
    % Merge data and simulation to restrict to plants common to both
    plant_table = innerjoin(plant_it,model_it,'Keys',{'id_plant'});
    plant_table = plant_table( ~any(isnan(table2array(plant_table)),2), : );
   
    % Call scatter plot and save to file
    filename = [ filestub sprintf('_%02.0f',t) ];
    titleString = [ titlestub sprintf(': period %2.0f',t) ];
    plotFitScatterEmissions( plant_table.Eit_model, plant_table.Eit_data, ...
        filename, titleString, inlogs, optfig, Elim )
end

end